Gravitational cooling of self-gravitating Bose-Condensates 
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ABSTRACT 

Equilibrium configurations for a self-gravitating scalar field with self-interaction are con- 
structed. The corresponding Schrodinger-Poisson (SP) system is solved using finite differences 
assuming spherical symmetry. It is shown that equilibrium configurations of the SP system are 
late-time attractor solutions for initially quite arbitrary density profiles, which relax and virialize 
through the emission of scalar field bursts; a process dubbed gravitational cooling. Among other 
potential applications, these results indicate that scalar field dark matter models (in its different 
flavors) tolerate the introduction of a self-interaction term in the SP equations. This study can 
be useful in exploring models in which dark matter in galaxies is not point-like. 

Subject headings: galaxies: evolution - galaxies; formation - galaxies: halos - dark matter 



1. Introduction 

One of the most fundamental problems in mod- 
ern Cosmology is the nature of dark matter. This 
problem itself has provoked an intense theoreti- 
cal search for realistic candidates able to drive 
the formation of the cosmological structure we ob- 
serve today (Sahni 2004). Among others, there 
are some candidates that have for long been con- 
sidered as exotic; a paradigmatic example is the 
case of scalar fields, which are the simplest parti- 
cles we can think of (Matos & Ureiia-Lopez 2004; 
Sahni & Wang 2000). However, no fundamental 
scalar field has ever been detected. 

Nevertheless, theoretical research on scalar 
candidates is a very interesting field in Cosmol- 
ogy nowadays, and many ideas have been pro- 
posed. One principal idea is that scalar fields 
can form Bose condensates, and then are able 
to form structure (Hu et al.2000; Arbey et al.; 
Matos & Urena-Lopez 2001). The quest for the 
kind of objects that are formed by scalar fields 
has been investigated in detail using numerical 



tools. The relativistic regime has been studied 
in (Seidel & Suen 1994; Balakrishna et al. 1998; 
Alcubierre et al. 2002; Alcubierre et al. 2003), 
for both complex and real scalar fields. 

Scalar fields are an example in which dark mat- 
ter is not point-like, hence no conclusions can be 
extracted from the widely studied N-body simu- 
lations for dark matter particles. Thus, we need 
to understand their properties in order to decide 
whether such properties could be useful to explain 
what we see in galaxies. On of the most attractive 
features of the scalar field model is that the con- 
struction of equilibrium configurations demands a 
smooth density profile at the origin, which is an 
advantage over the point-like dark matter particle 
models which have shown to be cuspy. 

In general, complex scalar fields can form sta- 
ble equilibrium configurations called boson stars, 
that are globally regular and whose energy den- 
sity is time- independent, for a recent review see 
(Schunck & Mielke 2003). Real scalar fields also 
have equilibrium configurations, that were dis- 
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covered in (Seidel & Suen 1991), and arc called 
oscillatons. The latter are also globally regular, 
but are fully time-dependent. As for their stabil- 
ity, it seems to be quite robust as far as numer- 
ical evolutions is concerned (Seidel & Suen 1991; 
Alcubierre et al. 2002; Alcubierre et al. 2003), 
but it may be that they are only long-lived 
(Page 2004). 

An important part of any scalar field model of 
dark matter is the scalar potential V {(()), which 
encodes in itself the self-interactions of the scalar 
field other than gravitational. The most popular 
scalar potential, and also the most studied in Cos- 
mology, is the quadratic one, V{^) = (m^/2)(^^ 
(or y(|(^|) = TO^|(/)p if the scalar field is complex), 
where the parameter m is identified as the mass 
of the boson particle. 

On the other hand, the inverse of the mass, 
i.e. the Compton length of the boson Ac = 
(in units in which ft = c = 1), becomes the 
natural unit of distance of scalar configurations. 
Actually, the most compact scalar objects have 
a radius of a few Compton lengths. In addi- 
tion, the boson mass also establishes the mass 
scale of the compact object, which is of the order 
~ mp;/m, where mpi is the Planck mass. Thus, 
a light scalar mass may provide of very massive 
objects(Seidel & Suen 1991; Harrison et al. 2002; 
Alcubierre et al. 2003; Schunck & Mielke 2003; 
Guzman & Ureha-Lopez 2004). However, scalar 
potentials more complicated than the quadratic 
one can provide of extra parameters which can 
add new features to oscillatons and boson stars. 

One simple instance of such more general scalar 
potential is that containing an extra quartic cou- 
pling in the form V{(p) = {fn^ /2)4>'^ + Ai„t?;>*/4, 
see (Colpi et al. 1986; Balakrishna et al. 1998; 
Arbey et al.) and references therein. In this case 
we speak of a self-interacting scalar field^ . 

The influence a quartic self-interaction can have 
in the dynamics of a cosmological complex scalar 
field has been widely studied in the literature, for 
recent developments see(Schunck & Mielke 2003; 
Arbey et al.). However, there are also some 
changes with respect to the free case. The self- 
interaction term can become important in situ- 



Thc quadratic potential is also called free potential, as the 
scalar field represents in this case a free relativistic particle 
of mass m with total energy E'^ = + m? . 



ations in which the scalar field has large ampli- 
tude, which may happen typically in the early 
universe, and in the interior of compact scalar ob- 
jects. For the purposes of this paper, it should be 
mentioned that a quartic term makes, for a given 
scalar field strength, a boson star larger and more 
massive than in the free ease (Colpi et al. 1986; 
Balakrishna et al. 1998; Guzman 2006). 

On the other hand oscillatons with self-interacting 
terms in the scalar potential have not been studied 
yet. The reason for this may be an intrinsic dif- 
ficulty in dealing with time-dependent relativistic 
configurations. 

Nevertheless, a curious point is that the New- 
tonian regime is the same for boson stars and os- 
cillatons. The relativistic EKG equations in this 
regime are equivalent to the so-called Schrodinger- 
Poisson (SP) system, see (Seidel & Suen 1994; 
Harrison et al. 2002; Guzman & Urefia-Lopez 2003; 
Guzman & Urena-Lopez 2004) and references therein. 

For a free scalar potential, the Newtonian equi- 
librium configurations can be as large as a galac- 
tic halo if the mass of the boson is very light 
(Arbey et al.; Guzman & Ureha-Lopez 2003). That 
the Newtonian regime is the astrophysically inter- 
esting one can be seen if one calculates the Comp- 
ton length associated to a (very light) boson mass 
of order 10~^^ eV; this is Ac ~ 1 pc(Arbey et al.; 
Sahni & Wang 2000; Matos & Urefia-Lopez 2001; 
Hu et al.2000). From the point of view of parti- 
cle physics, this boson mass is extremely small. 
Therefore, such a light boson will form relativistic 
objects which are small compared to typical galac- 
tic scales. However, in the Newtonian regime, 
gravity is weaker and allows the formation of much 
larger scalar objects. 

Within the Newtonian regime, there are also 
equilibrium solutions for different number of nodes 
of the wave function; solutions with nodes are 
called excited Newtonian boson stars and are clas- 
sified according to the number of nodes the wave 
function has. Zero-node states are called ground 
states, and configurations with nodes are generi- 
cally called excited states. Such excited states were 
also considered in the past as candidates for dark 
matter halos in galaxies (Lee & Koh 1996). They 
have the nice property of providing nearly flat ro- 
tation curves if the number of nodes is sufficiently 
large, see the example shown in (Arbey et al.). 
Unfortunately, those excited boson stars are 
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not stable, and they decay into a ground state 
in a time much smaller than the actual age 
of the universe (Guzman fc Ureiia-Lopez 2003; 
Guzman & Ureiia-Lopez 2004); the attractor be- 
havior of the ground state configurations is shown 
below in this paper for initially excited state equi- 
librium configurations. This fact rules out the 
possibility of excited boson stars as dark matter 
halos; nevertheless, they could still play the role 
of transient states in the evolution of scalar field 
configurations. 

It is interesting to note that the SP sys- 
tem has been widely studied in other fic;lds 
out of cosmology, see the interesting works in 
(Choi 1998; Choi 2002; Giovanazzi et al 2001). 
Also, they have been proposed to ameliorate the 
quasar's redshift problem in (Svidzinsky 2004). 

In this paper we want to study the case of scalar 
configurations with a self-interaction term in the 
Newtonian limit, that is, we shall study the SP 
system including a quartic self-interaction term. 
In doing this, we will learn about the properties 
that both a boson star and an oscillaton have in 
the weak field limit. 

Scalar field configurations with a quartic self- 
interaction have been proposed before to be 
the dark matter in galaxies, see (Arbey et al.). 
But, it is necessary to understand the dynam- 
ical properties of such configurations. In this 
sense, this paper extends the research already 
done in(Guzman & Urena-Lopez 2004) for the 
(free field) SP system. 

Assuming spherical symmetry, the Schrodinger- 
Poisson (SP) system of equations to be solved 
reads 

idri^ = -^a^(a;^) + [/V + A|V|'^ (1) 
dlixU) = xi^r- (2) 

We have used the dimensionless variables r — mt, 
X = mr, where t and r are the physical time and 
radial coordinates, respectively, and m is the mass 
of the boson. 

The wave function ^ is coupled to its own 
gravitational potential U, and the constant A 
represents the s-wave scattering length in the 
Gross-Pitaevskii approximation for Bose conden- 
sates (Gross 1963; Pitaevskii 1958). 

It is known that the SP set of equations with- 
out self-interaction has a fancy scaling property 



(Guzman & Ureiia-Lopez 2004) that simplifies 
the numerical research. It is possible to establish 
a similar scaling property when a self-interaction 
term is included. That is, Eqs. (1) are invariant 
under the transformation 

{r, X, U, ip, A} ^ {A-V, X-^x, X^U, X^ip, X'^A} 

(3) 

with A an arbitrary parameter. Unfortunately, 
there is a main difference with respect to the non- 
interacting case. In the latter, the solution of the 
SP assuming |'0(O)| = 1 generates the whole and 
unique branch of possible equilibrium solutions. 

For the new SP equations (1), it is necessary 
to find first all possible solutions that belong to a 
single value of the self-interaction coefficient A. If 
we now scale the latter using a certain A, we can 
generate a full new branch of solutions for A out of 
the one belonging to the old A. Though not as easy 
as in the free case, the scaling transformation (3) 
will still help us to understand the evolution of a 
scalar field configuration under an arbitrary value 
of A. 

There is also a complete set of physical quanti- 
ties that are well defined in the Newtonian limit. 
They obey the following scaling transformation 

{p, N, K, W, 1} {A^p, XN, X^k, X^W, X^i} . 

(4) 

where p is the density of probability, N is the in- 
tegral of p over the whole space, K and W are 
the expectation values of the kinetic and gravi- 
tational energies, and / is the expectation value 
of the self-interaction energy. The above quan- 
tities will help us to follow the evolution of an 
arbitrary scalar field configuration and establish 
criteria about virialization and equilibrium, 

This paper is organized as follows. In Sec. 2, 
we present the solution to the initial value prob- 
lem that will be the equilibrium configurations of 
the SP system. In Sec. 3, we describe the numer- 
ical methods and the boundary conditions used 
for the time evolution of the SP system. After 
that, in Sec. 4, we describe what we mean by 
virialization of non-equilibrium configurations and 
the late-time attractor behavior of the sequence of 
equilibrium configurations. Finally, in Sec. 5, we 
draw some comments and conclusions. 
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2. Equilibrium configurations 

For equilibrium configurations, we assume 
tiiat tiie wave function is of the form ip{T,x) = 
e*"^(/)(x), wiiere w is a free parameter. Ttiis as- 
sumption implies that the density of probabil- 
ity p and the gravitational potential U are time- 
independent, whereas the wave function evolves 
harmonically in time. Eqs. (1) and (2) become 

= 2x(?7-w)-F2A|(/)|20, (5) 
dlixU) = xcj>\ (6) 

The above system has to be solved under the con- 
dition of regularity at the origin ^(0) = dx(j)(Q) 
0. and isolation 0(x oo) — 0; we will also de- 
mand that U (oo) = 0. These boundary conditions 
determine in a unique manner the values of U (0) 
and Lu, which are the only free parameters of the 
solution. 

We made a numerical routine that solves 
Eqs. (1) and (2) using a shooting procedure that 
tunes the values of uj and U{Q) for a given central 
value (/>(0). In Fig. 1, we show the profiles of the 
initial wave function for some positive and neg- 
ative values of the self-interaction coefficient A. 
Also shown in there are the branches of equi- 
librium configurations for four different values 
of A. As expected, for A > there is nothing 
like a maximum indicating an unstable branch, 
which is always present (and typical) in rela- 
tivistic boson configurations (Colpi et al. 1986; 
Balakrishna et al. 1998; Guzman 2004). 

3. Numerical evolution 

In order to numerically evolve the SP system, 
we make a discretization of space and time, and 
approximate all derivatives using second order 
accurate finite differencing. The SP system is 
evolved one time step Ar using Eq. (1) to obtain 
a new wave function, and then we solve Eq. (2) to 
find the corresponding (new) gravitational poten- 
tial. 

We use an explicit time integrator to solve 
Eq. (1), as opposed to the fidly implicit method 
used and described in(Guzm;m & Ureha-Lopez 2004), 
where the problem of evolution was reduced to a 
linear system of equations. In the present case, 
due to the non-linear term in the Hamiltonian 



of Eq. (1), the reduction to a linear system of 
equations is not that simple. 

In any case, independently of the numeri- 
cal method used to integrate in time, we de- 
mand the evolution operator to preserve the num- 
ber of particles for an equilibrium configuration 
N = ^ \ip\x'^ dx"^ . With a modified version of 
the iterative Crank-Nicholson evolution method 
(Teukolsky 2000), we were able to reproduce the 
results found with the implicit method for the 
linear case, and confirm that the evolution was 
mass-preserving for a wide range of values of A. 

Boundary conditions. At every time step, 
Eq. (2) is integrated inwards, and thus we ap- 
plied the following boundary condition at the two 
outermost points of the numerical grid 

U{Xn-l) = -N{Xn-l)/Xn-l, (7) 
U{Xn) = -N(Xn)/Xn. (8) 

where n labels the outermost point of the numer- 
ical domain. The whole profile of U{x) is then 
found using a sixth order accurate Numerov algo- 
rithm (Koonin & Meredith 1990). 

As it is discussed in (Guzman & Urefia-Lopez 2004), 
the boundary conditions (7) and (8) are equiv- 
alent to impose the condition \ip(T,Xn)\ 
on the wave function; hence, we are forcing 
the system to remain in the computational do- 
main. Because we want to evolve systems 
out of equilibrium and allow the flow of par- 
ticles out of the numerical domain, wc imple- 
mented a sponge over the outermost points of 
the grid, which consists in adding an imaginary 
potential Vj{x) to the Schrodingcr equation (see 
(Guzman & Ureiia-Lopez 2004; Israeh 1981)). The 
expression we use for the sponge profile is 

Vj = -^Vo {2 + tanh [{xj - Xc)/S] - tanh (xc/S)} , 

(9) 

which is a smooth version of a step function with 

amplitude Vq, centered in Xc and width 5. the 
minus sign warranties the decay of the number of 
particles at the outer parts of our integration do- 
main, that is, the imaginary potential behaves as 
a sink of particles. It is also worth noticing that 
no term related to the self-interaction appears in 

^In full units, the total mass is given by M = (mp^/m)N, 
where wipj is the Planck mass. 
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Initial profiles of (|)(x) N vs r^j 




Fig. 1. — (Left) Profiles (f>{x) of equilibrium solutions of the SP for different values of the self-interaction 
coefficient A. As expected, the larger the coefficient A the more massive the equilibrium configuration is. 
(Right) Sequences of equilibrium configurations for different values of A. Each point in the curves represents 
a solution of the initial value problem of a total number of particles N and 95% radius rgs (the radius inside 
which 95% of the total mass is contained in). In this plot it is manifest that the bigger the A the less 
compact a configuration is. The filled circles indicate two configurations we use as tests for our numerical 
implementation. 



this conservation equation since we are assuming 
A is real. 

Tests. The obstacles our code must sort out are: 
i) the evolution of equilibrium configurations, in 
which the wave function oscillates with the definite 
frequency u found in Sec. 2, whereas the density of 
probability and the gravitational potential should 
remain timc-indcpendcnt; ii) the convergence of 
physical properties of the system (4). 

These two tests are shown at once in Fig. 2. 
On the left hand side, we show the convergence 
of the central density and the independence of the 
evolution on the location of the numerical bound- 
ary. We know that in the continuum limit, these 
configurations should evolve keeping p{t, 0) = 1 
for all times. It is shown that our approximation 
through finite differences converges to that value 
in the continuum limit. 

The central density shows oscillations whose 
amplitude converge to zero as the numerical grid 
is refined. We relate these oscillations to the re- 
sponse of the system to the intrinsic perturbation 
induced by the truncation error of the finite dif- 
ferencing approximation^. On the right hand side, 

^In (Guzman & Urcna-Lopcz 2004), we compared this type 
of oscillations for A = with a hnear perturbation analysis 



we show how the evolution operator preserves the 
number of particles all along the simulation. In a 
few words, what we show in Figure 2 is how the 
evolution method, the finite differencing and the 
boundary conditions work properly all together. 

Besides the expected results found for the cases 
marked with filled circles in Fig. 1, we also tried 
to evolve configurations located to the left of the 
maximum for A = —0.2. Such configurations are 
very compact, and their evolution resulted in a 
quick collapse and in the divergence of the cen- 
tral density p(r, 0). In principle, this would be 
the expected behavior for an unstable equilibrium 
configuration in the relativistic regime. However, 
we do not account with a clear definition of a col- 
lapsed object in terms of trapped surfaces of hori- 
zon formation during the evolution, so that it is 
not possible to draw accurate statements within 
the non-relativistic SP formalism for A < 0. One 
possibility would be to continue the evolution of 
the system with a relativistic code after certain 

and found excellent coincidence with the Fourier analysis 
of the numerical evolutions. In the present case, we deal 
with a non-lineaj: Schrodinger equation, and such analysis 
would say little about the cases with A 0. Nevertheless, 
we already have evidence that the truncation error could 
be responsible for this very regular oscillations. 
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Fig. 2. — Evolution of equilibrium configurations with different values of A. Figures on the left. For each 
case, the central value of p(t, 0) is shown for two resolutions: Aa; = 0.08,0.16 with the boundary located 
a.t XB = 50, and a third case with the boundaries located at Xb = 100 and resolution A = 0.08. The curve 
corresponding to the coarse resolution has been scaled as (/9o.i6(0) — l)/4+ 1, and the fact that it lies upon 
the value found for po.os indicates second order convergence of the solutions. Hence, we are running the 
simulations in the convergence regime of our discretization. The fact that the density found for Ax = 0.08 
with the boundaries at two different locations lies upon each other, indicates that the results are independent 
of the location of the boundary. Figures on the right We show the evolution of the number of particles, and 
the value of the virialization quantity 2K + 3/ + W, which in the continuum case should be zero. We 
see that its values are four times smaller if the spatial resolution is doubled. This clearly indicates second 
order convergence of this quantity to zero, which in turn means that we can recover what is expected in the 
continuum limit. 
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compactness is achieved, however such analysis is 
beyond the scope of this paper. 

4. Gravitational cooling and virialization 

The viriahzation and stabihzation of physi- 
cal configurations evolving according to Eqs. (1) 
and (2) with A = was first proposed in(Seidel & Suen ; 
and revised in detail in (Guzman & Ureha- Lopez 2004). 
The mechanism for relaxation consists in the ejec- 
tion of scalar field bursts, and the tendency of the 
system to settle down onto equilibrium configu- 
rations. This process was dubbed gravitational 
cooling. 

The criterion to determine whether or not a sys- 
tem was near an equilibrium state involves moni- 
toring the virial theorem relation 

2K + 31 + W = 0, (10) 

where K, I, W are the expectation values of the ki- 
netic, self-interaction, and gravitational potential 
energies, respectively, for the SP system (see for 
instance (Wang 2001)). In Fig. 2, we also show the 
oscillatory behavior of the quantity 2K -|- 3/ -|- W, 
whose amplitude converges to zero in the contin- 
uum limit. 

In Fig. 3, we show the virialization process for 
two different values of A. The initial profile is 
of the form 'ip{0,x) — cos^(a:;)e~^ , with a a 
free parameter. The reason why we choose such 
initial profile is that it provides a high initial ki- 
netic energy, and the process of virialization (if 
any) will be evident while calculating Eq. (10); 
other initial profiles of garden variety have been 
evolved, however the overwarmed initial configu- 
rations shown here present in a better way the 
relaxation process. The case A = is as simple 
as those found in the past (Seidel & Suen 1994; 
Guzman & Urefia-Lopez 2003; Guzman & Urefia-Lopez 
because it is possible to choose any profile with the 
confidence that the configuration will approach a 
rescaled equilibrium configuration (as shown in 
(Guzman & Urefia-Lopez 2004)), but we want to 
be sure that the same happens for A 7^ 0. 

Our expectations were fulfilled because of the 
manifest results in Fig. 3. All of the studied cases 
show a relaxation process towards one of the equi- 
librium configurations shown in Fig. 1. With these 
results, we confirm that the equilibrium configu- 
rations in Fig. 1, which are solutions of Eqs. (5) 



and (6), are indeed late-time attractor solutions. 

To finish with, we show some numerical runs of 
excited states for A = and A = 0.2 in Fig. 4. As 
expected, all excited configurations we tried are 
intrinsically unstable, and eventually settle down 
onto an equilibrium configuration. The initial ex- 
jUted equilibrium configurations start at the far 
tot)-right in the plots, and after a short time the 
systems approach the branches of ground state 
configurations shown in Fig. 1; the plot corre- 
sponding to the case A = is a reproduction of 
Fig. 16 in (Guzman & Urefia-Lopez 2004) with 
the appropriate redefined radius. The time of de- 
cay is of the order lO^rn"^, which is a quite small 
time scale compared to the age of the universe, 
even for small values of the boson mass^. 

5. Conclusions 

We have found stationary solutions to the 
spherically symmetric SP system of equations 
with a self-interaction term of the Gross-Pitaevskii 
type, according to the mean field approximation 
of Bose condensates. We found that equilibrium 
configurations exist that are virialized and stable 
under the action of small perturbations. 

For initial profiles with a high kinetic energy, 
we have shown that there are bursts of scalar field 
involving a considerable amount of matter of the 
order of half the initial mass and for stationary 
excited initial configuration the emission of mass 
is even higher; after a while, the system relaxes 
and virializes around an equilibrium configura- 
tion. This shows that equilibrium configurations 
are late-time attractors. We have made numerical 
experiments with less exotic initial profiles, like 
gaussians, and found the same tendency to viri- 
alize and to accommodate around an equilibrium 
juration. 

hat this shows is that a rather arbitrary ini- 
tial fiuctuation of self-gravitating scalar field, in 
the Newtonian regime, always evolves towards 
a virialized configuration, which should be com- 
pared to the relativistic case, in which some 
particular initial configurations lead to the for- 
mation of black holes, see(Seidel & Suen 1991; 
Balakrishna et al. 1998; Alcubierre et al. 2003; 

*Even for an ultra-light boson mass of 10~^^ eV, the time 
decay is of order ~ 100m~^ ~ 10^ yr. 
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A=0 A=0 




2000 4000 6000 8000 10000 2000 4000 6000 8000 10000 

T T 

Fig. 3. — Top panel We show the evokition of the initial profile -0(0,2;) = cos^(.T)e^^ /(3.959) A = 0. 

The system starts with a high kinetic energy, and after a relaxation period the central density p(t, 0) 
stabilizes around one, whereas the virial relation 2K + M + W oscillates around zero. The assymptotic 
behaviors of the the number of particles TV and total energy arc shown on the right plot. The initial mass 
is Mi = 3.617337 {rn?pi/m), and during the process around 42% of it was ejected very quickly. Bottom panel 
This is the case for A = 0.2, with an initial wave function profile ^l){0,x) = cos^(x)e~^^/^^'°^^)^. As before, 
the system also relaxes and stabilizes around an equilibrium configuration. This time, the initial mass is 
Mi = 3.98 (mp;/m), and the ejected mass is ~39%. In both cases the final configurations are nearly those 
shown in Figure 2. 
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Fig. 4. — (Left) The evolution of initial equilibrium configurations in excited states with A = 0. Such 
configurations start up with a given mass and radius (top-right in the figure) and after a while they approach 
one of the equilibrium solutions in the ground state through the emission of bursts of scalar field particles; 
the late time tendency is the approach to one of the points in the solid curve that contains all the possible 
ground state equilibrium configurations. This is a typical case of gravitational cooling. (Right). The same 
attractor behavior for A = 0.2. In both plots the solid lines correspond to those shown in Figure 1. For all 
these systems the quantity 2K + ZI + W was oscillating around zero with an amplitude converging to zero. 



Guzman 2004). Hence, a Bose condensate made of 
ultralight scalar field particles, as those proposed 
to be the dark matter, tolerates the introduction 
of a self-interaction term in the mean field approx- 
imation. 

We expect the results presented all along this 
work could be useful for models of scalar dark mat- 
ter, as some papers have explored the idea of in- 
cluding self-interaction terms (Matos & Ureha-Lopez 
Arbey et al.; Sahni 2004). We also hope this work 
will help others to decide whether such objects can 
be found in the cosmos. 

This research is partly supported by the bi- 
lateral project DFG-CONACyT 444-113/16/3-1; 
CONACyT grants 32138-E, 34407-E and 42748; 
PROMEP grants UGTO-CA-3 and UMICH-PTC- 
121; CIC-UMSNH-4.9 and Concyteg 05-16-K117- 
032. The runs were carried out in the Ek-bek 
cluster of the "Laboratorio de Supercomputo As- 
trofisico (LASUMA)" at CINVESTAV-IPN. 

REFERENCES 

Alcubicrrc, M., Guzman, F. S., Matos, T., Ni'mez, 
D., Ureha-Lopez, L. A. and Wiederhold, P. 
2002, Class. Quantum Grav. 19, 5017. 



Alcubierrc, M., Bcccrril, R., Guzman, F. S., 
Matos, T., Ni'mcz, D. and Urciia-Lopez, L. A. 
2003, Class. Quantum Grav. 20, 2883. 

Arbey, A., Lesgourgues, J. and Salati P., 2001 
Phys. Rev. D64, 123528; ihid2QQ2, 65, 083514; 
ibid 2003, 68, 023511. 

Balakrishna, J., Seidel, E. and Suen, W-M. 1998, 
; Phys. Rev. D58, 104004 

Choi, D. L 1998, Ph D Thesis, The University of 
Texas at Austin. 

Choi, D. I. 2002, Phys. Rev. A66, 063609. 

Colpi, M., Shapiro, S. L. and Wasserman, L 1986, 
Phys. Rev. Lett.57, 2485. 

Giovanazzi, S. O'Dell, D. and Kurizki, G. quant- 
ph/0010045 

Gross, E. P. 1963, J. Math. Phys. 4, 195. 

Guzman, F. S. 2004, Phys. Rev. D70, 044033. 

Guzman, F. S. 2006, Phys. Rev. D73, 021501(R). 

Guzman, F. S. and Ureha-Lopez, L. A. 2003, 
Phys. Rev. D68, 024023. 



9 



Guzman, F. S. and Urciia-Lopez, L. A. 2004, 
Phys. Rev. D69, 124033. 

Harrison, R., Moroz, I., and Tod, K. P. 
2002, math-ph/0208045; ibid 2002, math- 
ph/0208046. 

Hu, W, Barkana, R and Gruzinov, A. 2000, 
Phys. Rev. Lett.85, 1158. 

Israeli, M and Orszag, S. A. 1981, J. Gomp. Phys. 
41, 115. 

Lee, J-W and Koh I-G, 1996, Phys. Rev. D53 
2236. 

Koonin, S. E. and Meredith, D. C. 1990, Com- 
putational Physics (Addison- Wesley Publishing 
Gompany). 

Matos, T. and Urciia-Lopez, L. A. 2001, 
Phys. Rev. D72, 063506. 

Matos, T. and Ureiia-Lopez, L. A. 2004, Int. J. 
Mod. Phys. D13, 2287. 

Page, D. 2004, Phys. Rev. D70, 023002. 

Pitaevskii, L. P. 1958, Zh. Eksp. Teor. Fiz. 34, 
1240. 

Sahni, V. and Wang, L. M. 2000, Phys. Rev. D62, 
103517. 

Sahni, V. 2004, astro-ph/0403324. 

Schunck, F. S. and Mielke, E. W. 2003, Glass. 
Quantum Grav. 20, R301. 

Seidel, E. and Suen, W-M. 1994, 
Phys. Rev. Lett.6, 1659. 

Seidel, E. and Suen, W-M. 1994, 
Phys. Rev. Lett.72, 2516. 

Svidzinsky, A. A. 2004, astro-ph/0006024. 

Tcukolsky, S. 2000, Phys. Rev. D61, 087501. 

Wang, X. Z. 2001, Phys. Rev. D64, 124009. 



This 2-column preprint was prepared with the AAS lAT^^X 
macros v5.2. 



10 



